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ABSTRACT 


A simulation of collisional and gravitational interaction 
in the early solar system generates planets '^^1000 km in diameter 
from an initial swarm of kilometer-sized planetesimals , such 
as might have resulted from gravitational instabilities in the 
solar nebula. The model treats collisions according to 
experimental and theoretical impact results (such as rebound, 
cratering, and catastrophic fragmentation) for a variety of 
materials whose parameters span plausible values for early 
solid objects. Ad hoc sticking mechanisms are avoided. The small 
planets form in mio** yr, during which time most of the mass 
of the system continues to reside in particles near the 
original size. The relative random velocities remain of the 
order of a kilometer-sized body's escape velocity, with random 
velocities of the largest objects somewhat depressed due to 
damping by the bulk of the material. The simulation is 
terminated when the largest objects' random motion is of 
smaller dimension than their collision cross-sections, so 
that the "particle- in-a-box" statistical methods of the model 
break down. The few 1000 km planets, in a swarm still 
dominated by kilometer-scale planetesimals , may act as 
"seeds" for the subsequent, gradual, accretional growth into full- 


t"3 


sized planets. 



I . BACKGROUND 


Currently fashionable models for the formation of the 
planets require collisional accretion of planetesimals , 

Earlier theories had suggested that planet-sized objects 
might have formed as a result of turbulent vorticity which 
concentrated solid material at certain locations in the 
solar nebula (cf. Kuiper, 1951a) or as a consequence of 
gravitational instabilities which would accompany the 
flattening of solid material into a disk (Kuiper, 1951b) . 

More recently, Goldreich and Ward (1973) have shown that such 
effects might have produced planetesimals only a few kilometers 
in size in the region of the solar system now occupied by the 
terrestrial planets , 

How, then, were planet-sized bodies built from these plane- 
tesimals? Safronov (1972) suggested that mutual collisions 
and accretion produced larger bodies which then swept up the 
smaller ones within their gravitational cross-sections and 
scattered other planetesimals onto orbits which permitted 
later accretion. Safronov’s analytic approach required a 
decoupling of the evolution of random relative velocities 
of particles (i.e. orbital eccentricities and inclination) 
from the evolution of their size distribution. Safronov 
obtained the result that, once equilibrium is reached, 
relative velocities are comparable to the escape velocity 
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associated with that size particle which dominates the 
population, specifically the largest bodies in a reasonable 
power- law distribution. Many investigators have assigned 
that this condition applied throughout the growth process. 
Under this assumption Safronov modelled the evolution of the 
size distribution which yielded planet-sized objects whose 
final stages of accretion involved high velocity collisions 
which produced such observed properties as the distribution 
of obliquities and rotation rates of planets. 

This picture of planet growth has been refined by 
considering possible high velocity components of the colliding 
population. Weidenschilling (1974) hypothesized that Jupiter 
grew earlier than the terrestrial planets because low. tem- 
perature would have led to early condensation of more, and 
possibly stickier, solid matter. He suggested that Jupiter 
may have played a role in promoting accretion of terrestrial 
planets by inducing relative velocities among planetesimals, 
thus increasing the accretionary flux on the growing bodies . 
Weidenschilling (1975) and Kaula and Bigeleisen (1975) have 
proposed models in which planetesimals scattered by close 
encounters with Jupiter have dif ferent ' effects on each of the 
early terrestrial bodies and account for important differences 
in observed physical properties. 

Another source of enhanced relative velocities 
would have been resonances between the orbits of planetesimals 
and Jupiter. The possible importance of such resonances was 
stressed by Safronov (1972, p. 89) and by Kuiper (1974) . 
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However, a meaningful model of planet growth near resonance 
cannot be constructed by simply increasing relative velocities 
in Safronov's analysis because the induced relative velocities 
are size dependent (Greenberg, 1978). In fact, if we are 
to incorporate such sophisticated, but essential, mechanisms 
as orbital resonances into our ideas about planet growth, we 
must first devise a model that includes the coupled evolution 
of the size and velocity distribution. 

We have undertaken to develop such a comprehensive numeri- 
cal model for the evolution of already- formed swarms of plane- 
tesimals into small planets. Our aim is to include a wider 
variety of physical processes (e.g. resonances) and more detailed 
treatment of certain ones (e.g. low-velocity impact phenomena) 
than has been attempted before. The model has been developed 
and extended from an earlier model of Chapman and Davis (1975) 
which was intended to treat high-velocity asteroid -f ragmentation 
processes. The new model has provision for treating the entire 
range of planetesimal velocities and for treating both erosion 
and accretion, in addition to the catastrophic fragmentation 
processes already part of the original model. 

A fundamental element of this project is the careful 
evaluation — both experimentally and theoretically — of 
the nature of low- and moderate-velocity collisional inter- 
actions among solid bodies. This approach contrasts with 
previous models. For instance, the numerical simulation by 
Isaacman and Sagan (1977) ignores collisional mechanics 
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and simply assumes 100% efficient accretion of any particle 
impacting sequentially introduced accretion nucleii/ with 
no possibility of other particles sticking together. Very 
few experiments have been conducted previously of fragmentation 
and cratering processes at velocities much less than the 
hypervelocities that exist in the solar system today. It is 
expensive and time-consuming to investigate the entire range 
of velocities, projectile/ target diameters, material, properties , 
etc. for the whole suite of relevant parameters (rebound 
velocities, fragmental and ejecta size and velocity distri- 
butions, etc.). But through judicious selection of several 
low-velocity experiments combined with interpolation based 
upon physical principles, we can gain a much better under- 
standing of the low- velocity interaction of planetesimals . 

II. THE ALGORITHM 

In this paper, we report on the development of a computer 
simulation of the collisional evolution which includes the 
simultaneous variation of velocity- and mass-distributions 
with time. In this respect the model represents a significant 
advance toward an accurate portrayal of early events in the 
solar system. However, many important phenomena have yet 
to be included, and many parameters, sucij as the strength 
of the relevant material, can only be estimated. Nevertheless 
our algorithm provides a basis for a higher-order approxi- 
mation of, the evolution, and for study of the relative 
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influence of various phenomena and parameters . 

The population of particles in the collisionally evolving 
swarm is represented by a size distribution N(D), the number 
of particles per log increment in diameter D and over a par- 
ticular range of semi-major axes, a to a+Aa. This function 
is represented niomerically by the numbers of objects in each 
of the increments of log D over the size range under study. 

Each of these "size bins" has a particular variable value of 
the eccentricity e and of the inclination i associated with 
typical particles of that size. Ideally, one would like 
to have a distribution of e’s and i's associated with each 
size bin, but this would increase the number of dimensions 
to an unmanageable degree, although ultimately it may be 
possible to include at least some parameterization of distri- 
butions in e and i . The e ’ s and i ' s provide a measure of 
the random relative velocities among particles and the i's 
determine the thickness of the disk of particles in each 
size bin. We consider the events in a series of short time 
steps. In each time step, the probable number of collisions 
between particles in each pair of size bins is computed by 
a simple "particle-in-a-box" estimate: 

Number of oollisions - ''rel ^ section) x A time 

volume of space 


where Ni. and'Wa are the number of particles in two size bins, 
v^gj^ is the mean relative velocity computed from e's and i's 
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for both bins , the cross-section represents the geometrical 
cross-section enhanced by gravitational focusing, and the 
volinne is determined from a, Aa, and the disk thickness. 

Given the masses and velocities involved in each collision 
we determine the outcome of the collision in terms of 
resulting size distribution of ejecta, debris, or fragments 
and resulting relative velocities . The population distribution 
is then adjusted according to the number and outcome of 
collisions of each of several types disciissed below. The 
e's and i’s are modified by averaging in the new relative 
velocities of those particles which come out of collisional 
events. Time steps are chosen so that all changes in size and 
velocity in any one step are small. 

Collision Outcomes 

In order to implement this program we require some model 
of the outcome of collisions as a function of mass ratio 
and relative velocity at impact. For rocky materials 
experimental data are sparse compared with the wide range of 
masses and velocities required in our model. For this 
reason our model inevitably involves considerable extrapolation, 
which must be reconsidered as more relevant data become 
available. For now our algorithm divides the results of 
impact into four general categories, discussed in turn below; 

(i) Elastic rebound, (ii) Rebound with cratering of both 
surfaces, (iii) Shattering of the smaller body and cratering 
of the larger one, (iv) Shattering of both bodies. 
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(i) Elastic rebound . Impact at sufficiently low velocity 
between cohesive elastic bodies will result in rebound without 
chipping or cratering of either surface. Upon impact the 
surface of each body is depressed at a velocity v about half 
the impact velocity. The stress is alleviated primarily by 
propagation of a compressional wave at sound velocity c so that 
the maximum strain is ^-v/c. The maximum stress using the 
elastic modulus c^p would be cpv which must be less than the 
crushing strength S to prevent local fracture. For reasonable 
rock values this requires v < 5 m/sec or impact velocity 
< 10 m/sec. Indeed, our experiments show a transition' between 
no observable chipping at less than 10 m/sec to significant 
chipping at greater than 20 m/sec. 

The rebound velocity for basalt spheres at such low 
velocities is about 85% of impact velocity (Hartmann, 1978) ; 
for non-rotating irregular shaped rocks where substantial 
energy of collision goes into changing the rotation, this 
coefficient of restitution, c^^, is often less than 50% 

(Hartmann, 1978) . 

In our program the upper limit impact velocity for 
elastic rebound, v^, and coefficient of restitution are 
variable impact parameters. If the impact velocity (the mean rela- 
tive approach velocity augmented by the mutual gravitational 
acceleration) is less than the upper limit, the bodies separate 
at a velocity governed by the coefficient of restitution. If this 
separation velocity is less than the mutual escape velocity , the 
particles combine to produce a new particle whose mass is 
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the sum of their masses and with e and i dictated by the in- 
plane and out-of-plane components of the mean velocity of 
the center of mass of the two colliding particles with respect 
to ' circular orbits . If the rebound velocity is great enough 
to permit mutual escape, the particles remain distinct with 
their e*s and i's changed to reflect the change in relative 
velocities. - 

For bodies with regolith surfaces or with an intrinsically 
weak nature, such as primitive carbonaceous chondritic material, 
the upper velocity limit, v , for regime (i) might be 

practically zero. This type of material might well be 
representative of early solar system rocky condensates so 
collisions in regime (i) might never have occurred. On the 
other hand the treatment of this regime is incorporated 
into the program to permit flexibility in the types of material 
which can be studied. Conceivably, depending on heating 
mechanisms, early material might have achieved characteristics 
of hard rock shortly after accretion. We know that rocky and 
metallic bodies exist today. 

(ii) Rebound with' cratering of both bodies . If the rebound 
limit of (i) is exceeded but neither body is catastrophically 
disrupted, the surfaces of the impacting bodies will be 
locally damaged, e.g. chipped or cratered. (Hereinafter we 
shall refer to all such local damage as "cratering".) We 
know of no experimental results which give data specifically 
for this regime. On the other hand, there exist experimental 
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results in which high-velocity projectiles deliver kinetic 
energy to the surfaces of semi-infinite targets , thus producing 
craters. The mass excavated by these events can be approximated 
as a constant, K, times the kinetic energy delivered. Marcus 
(1969) summarizes results of Gault which give this constant K 
as AjS X 10 cgs for basalt, 'V'1.5 x 10“^ cgs for "weakly 
bonded quartz sand", and 'v2 x 10~® for sand. In our. algorithm 
we assume that half the kinetic energy of impact is delivered 
to each body in the colliding pair and we input some value 
of K within the plausible range to evaluate the mass cratered 
from each. 

The cratering process removes mass from each body, 
whereupon the bodies rebound at some fraction, of their 

impact velocity. This modified coefficient of restitution is 
less than the coefficient used in case (i) due to the' loss of 
energy in cratering and' mass ejection at the impact site. 

Indeed, in our experiments we have observed cases where a 
basalt ball or irregular igneous rock is fired into a rock 
target at 19 to 26 m/sec and undergoes minor cratering with 
•modified rebound; in the two measured cases with basalt balls, 
the cratered projectile rebounded from impacts at 26 and 29 
m/sec with velocity 0.73 and 0.90 times the normal rebound 
velocity, respectively. For weak materials the modified 
coefficient of restitution might be as low as 0.001. Just 
as, for case (i) , the computer program checks whether the 
•rebound velocity permits separation or accretion of these 
bodies . 
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What can we say about the velocities of the crater ejecta? 
The best available data on crater ejecta velocity was given 
by Gault/ Shoemaker and Moore (1963) . They provided a plot 
estimating the cumulative mass of ejecta vs. velocity for 
cratering into basalt at velocity '^'6 km/sec (Fig. 1, heavy 
line) . This result must be extrapolated over a great range 
of velocities/ materials, and mass ratios. Therefore, we 
have introduced a simplified version of their result (Fig. 1, 
dashed line) . Features such as the high velocity ejecta jet 
have been neglected as being too detailed for our degree of 
extrapolation. Marcus (1969) applied the same simplified 
curve for basalt to impacts in sand or regolith. His estimate 
was unreasonable because combined with the large value of K 
for sand, it gave the ejecta more kinetic energy than v/as put 
into the system by the impact! In fact, Stoffler, ^ (1975) 

show that ejecta from craters in sand at 6 km/sec travels 
'\zlO~'* times as far as for a comparable event in basalt. Hence 
we hypothesize that the velocities for sand can be represented 
by shifting the curve for basalt leftward by a factor of 10“^ 
as shown in Fig. 1. A major problem is applying these crater 
ejecta data to craters formed at lower velocities. The 
curve may be relatively independent of impact velocity if the 
percent of impact energy going into ejecta kinetic energy does 
not vary with velocity. Applying this argument, Gault 
et al. (19,63) estimated that impacts into basalt at tens of 
km/sec would give curves close to their 6 km/sec result. 
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We have performed experiments at much lower impact velocities; 

6 m/sec impacts in vacuum into fine rock powders give ejecta 
velocities consistent with the estimate for sand in Fig. 1. 

For real early solar system materials the relevant curve 
for cumulative fraction (f) of ejecta with velocity greater 
than v may resemble the intermediate curve: 

f % c v^^/** where c„-: = 3 x 10® cgs. The coefficient 
ej 

c . is an input parameter for the algorithm. The fraction of 
ejecta escaping from the parent body is given by the value 
of f which corresponds to the parent ’ s escape velocity . 

If that value of f is less than 1, the bulk of the escaping 
ejecta barely escapes so we take the e and i values to be 
the same as those of the parent body. Here by "parent body" 
we mean the body from which the pieces were ejected or, if 
they accrete one another, we mean the combined body. If 
the escaping fraction is unity, the ejection velocity before 
escape is taken to be the value at which the f vs. v function 
intercepts f = 1. 

How are the ejecta distributed in size? Based on experi- 
ments, observations of natural fragments, and earlier literature, 
Hartmann (1969) found a size distribution N m“^/® where 
N is the number of particles with mass greater than m. This 
-2/3 power law applies to cases in which- the locally damaged 
region receives barely enough energy for breakage and 
ejection. The largest piece has a mass given by setting 
N = 1. The constant must have a value (M/2) , where M 
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is the total mass of. escaped ejecta, in order to conserve 
mass. For cases in which an excessive energy density is 
applied to the damaged region, Hartmann found that the power 
law exponent is closer to -1. This result pertains for 
cases of impact velocity greater than the speed of sound in the 
material (i.e. "hypervelocity" impact), because energy 
propagates away from the impact site slowly compared to 
the rate of impact energy delivery. 

(iii) Smaller body shattered and larger one cratered . 

If sufficient energy is imparted to an entire body, 

it will fragment catastrophically, rather than experience 

merely local cratering. What are the criteria for catastrophic 

fragmentation? Most collision experimentation has dealt 

with hypervelocity cratering in semi-infinite 

targets. In such cases, the target is damaged only locally 

while the "bullet" undergoes super-catastrophic failure. 

Only a few experiments have been performed with targets 
small enough to yield results near the catastrophic limit. 
Theoretical evaluation of impact strength (Harris, 1975) 
has been unsuccessful because the processes involved are so 
complex (superposition of surface-reflected seismic waves, 
energy loss to heat and rotation, etc.). We have performed 
experiments with both finite and semi-infinite targets at a 
range of velocities (3-300 m/sec) . Catastrophic failure 
occurs if a critical energy per unit volume is delivered 
to a body by an impact. In most collision experiments there 
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is a sharp transition over a, narrow range of energy densities 
from minute local cratering to massive body fracture. For 
rocky materials and water ice these values are estimated to 
be 3 X 10^ and 2 x 10® ergs/cm^ , respectively (Hartmann, 1978; 
Greenberg e^ / 1977). For dirt clods, which may have cohesive 
strength comparable to early solar system solids, Hartmann (1978) 
finds a value of 'x^lO®. This parameter ("impact strength")' has 
the dimensions of strength (supportable force/area) but is 
conceptually distinct. Experiments by Moore and Gault (1965) 
and by Fujiwara ^ (1977) confirm the impact strength 

for basalt even at much higher impact velocities (1-3 km/sec) . 
These results indicate that impact strength is independent of 
velocity. 

Further experiments are needed to show how energy is 
partitioned between colliding objects, and whether impact 
strength depends on object size. In our model we assume 
tentatively that half the kinetic energy of impact is 
delivered into each body, and that strength is independent 
of size. Hence, with increasing energy of collision, the 
smaller body of a colliding pair will shatter before the 
larger one. Planetary cratering is generally in this 
category. Studies of such impact events usually emphasize 
the cratering process and ignore the fate of the shattered 
projectile. In our study, where collisions between bodies 
of comparable masses must be considered as well as between 
bodies of very different masses, we must consider the debris 
from the smaller body as well as the crater ejecta. 
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How are the pieces of the shattered body distributed? 


Hartmann (1969) noted that, just as crater ejecta, such 
debris follows a cumulative size distribution ‘of the form 


N = Cm where C and b are constants; b varies from a 
value of about 2/3 for cases where fragmentation energy is 
minimal to about 1 where large amounts of excess energy are 
delivered. In order to construct a computer algorithm we 
■needed to quantify Hartmann's observation. We may estimate the 
mass of the largest fragment by taking N = 1 which yields 


_ ^1/h 


m 

max 
ments as M = 


Integration gives the total mass of the frag- 


■(«) 


1-b , 

%ax ' ^ equals the mass of the 


shattered body . - We may solve for b and C in terms of m and 

max 

find b = { + 1}-' and C = All we need in 

order to determine the size distribution is a way to calculate 

"'max* Note that if varies from M/2 (barely catastrophic) 

-towards 0 (super-catastrophic) , b varies from 2/3 to 1 in 

agreement with Hartmann's observation. Fujiwara ^ ar. (1977) give 

"'max^^'^ "" ^ (Energy/mass)-^ •2'* for basalt; for other 

materials we scale the coefficient so that m = M/2 at 

max 

the critical energy density. 


The velocity of the fragments with respect to 
the impact site is computed such that all have the same 
speed and their kinetic energy is 50% of the energy 
delivered to the shattered body by the impact. This value is 
consistent with our experimental results for basalt and other 
igneous rocks at. impact velocities of tens of meters per second, 
although we are neglecting any 'high velocity "tail" in the 
distribution. (In reality there must be some distribution of 
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debris velocities, but no relevant experimental data yet 
exist.) The algorithm checks whether the debris' velocity 
is sufficient to escape the larger body ' s gravitational 
field. If it is not, they fall back and accrete; otherwise, 
they escape and add to the numbers of smaller particles 
with their e's and i's averaged in with those of the 
pre-existing small particles. 

(iv) Both bodies shattered . If the energy of collision 
is sufficiently great, both bodies will shatter. Again, 
we assume that half the energy goes into each body. The 
fragment size and velocity distribution for debris from each 
body is computed as for the smaller body in (iii) . The total 
kinetic energy is compared with the grayitatjLonal 'binding' energy . 
If sufficiently small, the debris fall back together; 
otherwise many small particles are created. 

Fig. 2 summarizes our treatment of the collision outcomes 
as a function of mass ratio and impact velocity. Fig. 3 
illustrates the change in particle mass for each outcome 
category. 

Re-distribution of Sizes 

The^ size distribution can be changed in two distinct ways 
in any time' step. First, the number of particles in each 
size bin may be changed. For example, catastrophic fragmentation 
of a large body removes one body from its size bin and adds 
many bodies "to smaller size bins according to the pov/er law 
distribution. .The second mode of re-distribution is 
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more subtle. Bodies in a given size bin may change mass, as 
in accretion or cratering erosion, by increments too small 
to effect a transfer into a new size bin. Such incremental 
changes are crucial to an evolution model, especially if 
accretional growth is expected. This requirement contrasts 
with the asteroid collision model of Chapman and Davis (1975) 
in which the dom.inant collisional process was assumed to 
be catastrophic fragmentation. In order to account for 
incremental changes in mass we adopt the following procedure. 

The number of particles in each size bin is assumed to be dis- 
tributed uniformly in log D. During a time step the average 
change in mass is computed for particles in each bin. This 
shift in mass moves those particles at one end of the size 
bin into the next bin. The mass of particles shifted into 
the next bin is calculated and is used to compute a change 
in the number of particles in that next bin, in such a manner as 
to conserve mass. The e and i characteristic of the former 
bin are averaged into the e and i of the new bin. 

Growth of the bodies in the largest size bin by accretion . 
might shift some small portion of the bodies into a new 
largest size bin in each time step. In general, our algorithm 
suppresses this transfer until a sufficient mass change accumu- 
lates that the number of bodies transferred into the new bins 
is greater than one. Otherwise meaningless infinitesimal 
fractions of bodies would be placed into the large size bins. 
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Actually, the choice. of unity for the lower limit on numbers 
in the new bin is arbitrary. A fractional number of 
particles in a size bin can have meaning since it is 
the number per size increment per semi-major axis increment. 

(Both of these increments are arbitrary with only the 
restrictions that size bins be sufficiently narrow that 
no important structure of the size distribution is 
neglected and that Aa be small enough that eccentricity 
and inclination give a good estimate of random velocities 
over the entire a range.) We have experimented with various 
values of the required number for opening a new largest bin 
and find that the characteristics of the evolution are 
relatively independent of this choice. 

A different criterion for populating a new largest size 
bin 'is needed when the largest bodies are so large, and 
accrete so efficiently, that their masses increase very 
fast compared with their numbers. As they accrete one 
another according to the formalism of particle- in- a-box 
statistics, their numbers may decrease faster than smaller 
bodies grow to replenish their numbers. Hence, they may 
grow by an amount greater than the increment between 
size bins before the criterion described above permits 
a transfer of bodies into the next bin. For this reason we 
permit the transfer if the mass change is a significant fraction 
of the bin width, even if less than the normally required number 
is transferred into the new largest bin. 
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When the distribution reaches a stage where there are 
only a few bodies per A.U. over a significant range of the 
largest size bins,- the system begins to be dominated by the 
statistics of small numbers which our program is not designed 
to treat. The final stages of planetary accretion are 
simply not amenable to the particle- in-a~box approach. 

However, as we shall show, our algorithm does work over an 
evolutionary period in which thousand- kilometer bodies 
are produced from a population originally all '^<1 km in 
diameter . 

Re-.distribution of Velocities 

The orbital eccentricity and inclination represent the 
in-plane and out-of -plane components of the random relative 
velocity of particles with respect to purely circular Keplerian 
orbits. Collisions are assumed to result dominantly from 
these' -random motions, rather than from differential Keplerian 
velocities bringing bodies within their collisional cross- 
sections. For collisions between bodies from different size 
bins, the mean relative approach velocities are computed from 

both sets of random velocities . The mean in-plane and out-of-plane 
velocity components of the center-of-mass of the two colliding 

bodies are also computed. The collision outcome is then deter- 
mined in the center-of-mass reference frame. Velocities 

of any debris, ejecta, or rebounding particles after escape 
are then added to the mean center-of-mass velocity. The 
mean in- and out-of-plane velocity components of escaping 
material are computed assuming isotropic escape with respect 
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to the center-of-inass of the colliding system. Any newly 
created particles are distributed into appropriate size 
bins and their velocities are root-mean- square averaged 
with the random velocities already associated with particles 
of their size, yielding a corresponding adjustment of e and i. 

In general collisions tend to damp the random velocities, 
although this is not always the case. For example, a slow 
moving body may hit and rebound from a fast moving large one. 
Even if the coefficient of restitution is significantly less 
than one, the small body may gain kinetic energy. Safronov 
(1972) pointed out another mechanism that tends to increase 
random velocities: the gravitational interactions of close 

approaches. With a number of approximations and assumptions 
he performed an analysis which indicated that an equilibrium 
between this stirring effect and collisional damping would 
yield random velocities on the order of the escape velocity - 
of the dominant size particle. 

In our program we numerically compute the change in 
random velocities in each time step due to gravitational 
stirring. Just as for each size particle we consider the 
probability and consequences of collisions with each other 
size particle, so we also consider the gravitational stirring 
as it passes through the field of particles of each other size. 
The ultimate source of gravitational stirring is the relative 
velocity between particles due to their differential Keplerian 
periods (cf . Safronov, 1972) . Gravitational interaction rotates 
the relative motion so that a non-circular orbit is generated.. 
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This randomizing of Keplerian shear is modelled in the following 
way. If a mass m^ moves past another mass m 2 , its velocity 
changes due to the gravitational interaction by an amount 
6v perpendicular to the initial relative velocity, v. This 
change is given by 


'S V = 


m2V 


mi + m2 


sin 2x 


( 1 ) 


where 


sin X 


= ( 1 . 


p2^4 


—h 


G (nil ^2 ) 


(2) 


P is the "impact parameter” {cf.. Ward 1976). If mi 
moves through a field of particles of mass m 2 , the mean 
square change in velocity per unit time is given for the 
planar case by 



(3) 


where a is the number of particles of mass m 2 per area of 
the disk. 

If the velocity, v, is due to Keplerian shear, v = 
where n is orbital mean motion. This substitution and 
integration from small to large P yields approximately 


dv^ 

dt 


qu/z^i/3 ^ 


m" 


a ^ 
9/J 


(4) 


We may show that this stirring model is consistent 
with Safronov's result by considering a simplified case with 
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all particles of equal size (mass m and diameter D) . From 
Ward (1976) we get for the time scale for velocity damping 
due to energy loss in collisions in our notation 


X = 


/l 4 

\y dt / TrnD''3 o 


(5) 


where 3 is the fractional energy loss per collision. Thus 



damping 


~ nD^gav 


( 6 ) 


From (4) for mj = m^ : 
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In equilibrium (. 6 ) and (7) are balanced so 
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where = /4Gm/D is the escape velocity and p = 6 m/( 7 rD^) 
is the material density. Taking p 'X' 3 gm/cm^ , 3 .a- 0.8 and 

n a. 2ir/yr we obtain v 'v 7 v . 

e • 

The dependence of equilibrium random velocity on n~V^ 
is worthy of note. Over .the entire solar system the 
random velocity varies by less than an order of magnitude. 
But beyond the distance of Pluto the random velocities would 
have been significantly higher than the particles' escape 
velocities, a factor which may have inhibited accretion at 
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such great distances, from the sun. Moreover, for material 
around -other stars with different masses, Kepler's third 
law gives a different mean motion for a given distance. 

The relation between mean motion- and random relative 
velocities would limit the region in which planets might 
form. If this region fails to overlap the region in which 
temperature and pressure permit condensation, planet formation 
may be prohibited! 

The analytic approach to velocity determination requires 
such assumptions as a simplified and constant mass distribution 
and velocity equilibrium. Our numerical approach requires none 
of these assumptions . We simply compute the change in velocity 
and mass distribution during each time step. Besides collisions, 
which primarily damp relative random velocities, and randomization 
of Keplerian differential velocities, we also take into account 
the rotation of random velocities due to gravitational encounters 
which can convert in-plane and out-of-plane motion from one 
to the other, thus partitioning random velocities between 
eccentricity and inclination. 

III. NUMERICAL EXPERIMENTS AND EESULTS 

Since knowledge of the relevant initial conditions, as 
well as material properties, is minimal, we must regard our 
computer . simulation as a means of testing for the range of 
planetesiraal conditions which lead to planet building. Does 
collisional evolution lead inevitably to growth of large 
bodies, or are very special initial velocities, mass distribu- 
tions, and materials required? We have begun -to test for the 
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generality of planet . growth by selecting some plausible 
starting parameters and in subsequent numerical experiments 
varying these parameters to the limit of their reasonable 
range. So far, indications are that 'vlOOO-km bodies grow 
from 1-km bodies for a wide range of parameters and initial 
conditions. And they grow fast, on time scales of a few 
tens of thousands of years or less . For evolution beyond 
this stage we find that the random motion is too small to 
justify the "particle-in-a-box" statistics. These and other 
results will be discussed after the details of the various 
numerical experiments are presented. 

Experiment 1 (Figs. 4, 5 and 6) : Nominal Parameters 

We have begun our numerical experiments by considering 
evolution near semi-major axis a = 4 x 10^^ cm (about 2.7 AU) 
over a range (Aa) of 8 x 10^^ cm (0.05 AU) . We take the interval 
between size bins to be a factor of two in diameter. Initially 
we assume all bodies to have diameter 1 km as suggested by 
Goldreich and Ward's (1973) gravitational instability calculations 
We take 10^^ as the initial number of such bodies, because, 
for a material density of p 3 gm/cm^ , this niimber gives a 
surface density of the particulate disk of about 8 gm/cm^ , 
the value used by Goldreich and Ward in their calculations 
(comparable to the surface density computed by "spreading 
out" the mass of the present terrestrial planets over their 
portion of the present solar system; see also Lecar and 
Franklin, 1973) . 
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For initial random velocities we selected values a few 

times the escape velocity of the 1 km bodies. A choice in this 

range seems appropriate in light of our discussion of equilibrium 

random velocities. Gravitational stirring would prohibit much 

lower velocities because v would be raised to several times 

v^ on a time scale “vlOOO years according to eqn. (7-) . A 

much greater initial velocity than we selected might lead 

to shattering and comminution of debris rather than planetary 

growth. (This occurred quite dramatically when we performed 

one run with weak material in which the initial velocity 

was about 20 v^.) If this occurred in nature, one of two 

outcomes might result; (a) The comminuted debris might be 

removed by solar wind pressure, inhibiting planet growth by 

lowering the available mass or (b),if the comminuted debris is 

not 'removed, the material would reaccumulate into 'v-1 km 

sized bodies by the Goldreich-Ward process. Because we know 

that planets ultimately did grow from that stage, at some 

point the velocities must not have been too much greater 

than a few times v . 

e 

* 

The initial eccentricity and inclination for the 
planetesimals were each taken to be 5 x 10"“* which corresponds 
to random velocities of about 700 cm/sec, about nine times 
the escape velocity,. We found that the .first stage in the 
evolution is predominantly characterized by damping of this 
velocity down to less than half the initial value even 
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before any accretional growth takes place. Planet growth 

occurs independent of the choice of any initial velocity between 

4 and 9 v . 

e 

In the first numerical experiment parameters were 
selected to approximate a. material somewhat more loosely 
bonded than basalt, but more cohesive than merely gravitationally 
bound sand. (Parameters used for all experiments are summarized 

in Table 1.) The choice of K and c . means that about 3% 
of the impact kinetic energy goes into ejecta kinetic energy. 

The selected impact strength, S = 3 x 10^ ergs/cm^ , is 
perhaps somewhat too high to be consistent with the idea of 
weak early solar system materials, but later tests showed 
it doesn't seem to affect the evolution in a crucial way. 
Moreover, weak, loosely-bonded surface material does not 
necessarily imply that S is proportionately low. Conceivably, 
the interior would be packed more densely giving substantially 
greater resistance to catastrophic fragmentation than the 
surface properties would indicate. Also, such intrinsically 
weak material might be ineffective at propagating seismic 
energy of impact through its volume, thus inhibiting dis- 
ruption. 

The resulting evolution of the size distribution is shown 
in Fig. 4. The size distribution is shown near each time 
step at which bodies are placed in a previously unoccupied 
size bin. At each time step we output a matrix of outcomes 
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of collisions between bodies in each pair of size bins- For 
experiment 1, the matrix remained nearly constant throughout 
the evolution and is schematically shown in Fig. 5. Evolution 
of the eccentricity (in-plane random velocity) distribution 
is plotted in Fig. 6. (The inclination distribution is generally 
quite similar, within a few percent.) 

For the first 10900 yrs, this evolution consists of 
collisions between 1 km bodies which crater one another but 
rebound and escape from one another. A small amount of crater 
ejecta escapes in each interaction, creating a distribution of 
small particles which in turn crater and escape one another 
and the 1 km bodies. The erosion of mass shifts some 1 km 
bodies into the 500m size bin, and these in turn erode into 
smaller sizes. In this manner, all the smaller bins are 
populated albeit with only a small fraction of the total mass. 
About 0.1% of the total system mass is lost beyond the smallest 
size bin (30 m) and is ignored in our program. The eccentricity' 
of 1 km bodies damps down to about 0.00023 ('^4 v^) before 
accretion begins. The eccentricities of the smaller bodies damp 
down much more slowly because the bulk of the mass is in 
t*odies much larger than themselves so gravitational stirring 
is more effective relative to collisional damping. This result 
may seem counter-intuitive to people who think of small bodies 
as generally being more susceptible to drag due to their large 
area/mass ratio. 

Once the relative velocities of the 1 km bodies become 
low enough, the bodies begin to accrete one another. The sub-km 
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bodies never slow down enough to be accreted by the 1 km bodies. In 
fact, some of them have their random velocities pumped up by gravi- 
tational interaction with larger bodies . Bodies in diameter range 

2 to '^100 km accrete everything smaller that hits them. Bodies 
greater than '^^100 km accrete any impacting bodies after shattering 

them. Their greater gravitational cross-section permits more 
efficient accretion and hence introduces the reverse curve 
slope to the size distribution. A 500 km body is produced 
about 10,000 yr after accretion begins. By this time the 
statistics of small numbers must be important so our particle- 
in-a-box algorithm becomes invalid. 

The most striking feature of this evolution is that most 
of the mass of the entire system remains in 1 - 4 km size bodies 
even after 100 - 500 km size bodies are produced. (This fact 
is made evident by comparison with a line in Fig. 4 whose 
slope is- such that 8 times as many particles are in each 
succeedihgly smaller box. A line of such slope represents equal mass 
per size bin.) The random velocities of the larger objects are 
damped by the 1 km bodies which appear to them as a dense, viscous 
medium. The random velocities of all bodies are quite low, on the 
order of the 1 km bodies ' escape velocity as we would 
expect, because these bodies dominate the population. 

We conclude that in about 20,000 yr, a disk of 1 km bodies 
evolves to include a small number ('^^25 per A.U.) of- 300 - 700 km 
bodies in their midst. Such a small number of large bodies 
might form the seeds for subsequent growth of a few 
planets . 
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Experiments 2 & 3 (Figs. 7 & 8) ; Suppressed Ejecta Velocities 
In order to demonstrate the minimal role played by crater 
ejecta in the evolution/ we show {Experiment 2) a run executed 
with a somewhat different algorithm for ejecta distribution 
which effectively prevents escape of ejecta from its parent 
body. All ejecta were assumed to depart the surface with a 
velocity of 0.005 times the impact velocity, intermediate 
between the hypervelocity results of Stoffler et al. (1975) 
and of Gault et (1963) . This rule, in effect, prohibits any 

ejecta from escaping its parent body. The resulting evolution 
is shown in Fig. 7. The results are virtually unchanged from 
the previous case for sizes 1 km or greater. The need to 
follow evolution over a smaller number of size bins permitted 
simulation over a greater time range than in Experiment 1 for 
the same given computer time limit. This resulted in creation 
and growth of 1000 km sized bodies after 24,000 yr, although 
the small-number statistics after the creation of 500 km objects 
give us little confidence in the validity of subsequent evolution. 

This evolution was also simulated with a somewhat narrower 
value for the size bin interval, a factor of 2 ^ 1 ^ instead of 2 
(Experiment 3) . Again, the results (Fig. 8) are practically 
unchanged from those of Experiment 2 (Fig. 7) , indicating that 
the size bin interval of a factor of 2 is sufficiently fine 
to model evolution adequately. 

Experiment 4 • (Fig. 9) ; Large Initial Population 

In order to explore the importance of our choice of 
surface density of material in the disk, we ran one simulation 
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with 100 times as many 1 km objects as in the first run. The 
characteristics of the evolution (Fig. 9) were virtually identical 
with previous runs with one striking exception; the time scale 
was contracted by a factor of 'x^lOO. Creation of small particles 

from crater ejecta began almost immediately. Velocities of 
1 km size bodies damped to '^-4 v^ in 109 yr, at 
which time accretion began. A 500 km body was produced at 
t = 157 yr and a 1000 km body at t = 179 yr. At each stage in 
the accretion/ (i .e. whenever a particle was created in a new 
bin) the eccentricity and inclination distribution was similar 
to that of the previous runs. 

Experiment 5 (Figs. 10,11, and 12); Solid Rock 

Our next experiment gives some indication of the importance 
of the assumed material properties. As an extreme case we 
assume, the material to have impact properties of solid, 
cohesive, competent rock. (See Table 1) . The coefficients 
of restitution c^^ and c^^ were both taken as 0.86, Hartmann's 
(1977) value for smooth basalt balls. For initial conditions, 
the same number of 1 km bodies and e and i were used 
as before. The results are shown in Figs. 10, 11 and 12. 
Initially, the 1 km bodies rebound after impact with one another. 
No cratering occurs. Their random velocity damps down for 
109,500 yr until e i 1.4 x 10“ S equivalent to random 
velocity of v^/4- At this point, the energy lost in a 
collision is sufficient to prevent escape after rebound. 
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so accretion occurs. {The newly formed body is assumed to 
acquire strength properties comparable to the original material 
before the next impact. Since this is probably unrealistic, 
this solid rock case is clearly an extreme.) In a matter 
of a few years, bodies as large as 64 km are formed.. While 
bodies between 1 and 32 km continue to rebound and accrete 
gravitationally, ‘impacts into 64 km bodies are accelerated to 
sufficiently high velocities that cratering occurs, producing 
ejecta. As before, the ejecta never constitutes an important 
fraction of the mass in the system. Ejecta velocities 
are so high that the particles subsequently rebound without 
accretion from any bodies smaller than 64 km. The 1 km 
and larger bodies continue to accrete one another. Once 
bodies greater than 64 km are formed, they accrete any 
smaller impacting objects after shattering them. Within 
a hundred years after first accretion, bodies hundreds of 
kilometers in diameter are produced. 

One striking property of this evolution is the hump in the' 

size distribution (Pig. 10) at about 200 km. This may be an 

artifact due to the fact that, for the low approach velocities 

developed in this evolution, the formal gravitational collisionai 

cross-section of a body greater than about 300 km exceeds its 

gravitational sphere of influence. The two body equations of 

motion used tp compute the gravitation cross-section 
are invalid outside the sphere of influence, where 

solar gravity dominates. To model this effect in our algorithm, 

the cross-section is cut-off at the sphere of influence. Hence, 
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there may be a discontinuity in the efficiency of collisions 
for bodies greater than o-300 km which would act as a dam 
slowing further growth and creating the 200 km hump. 

Actually, the very low eccentricities and inclinations’ 

can be shown to invalidate evolution for this solid rock model 

for any bodies greater than about 40 km. These e's and i's 

correspond to excursions in distance vea = 10^ cm from 

circular motion. The gravitational cross-sections for bodies 

larger than 40 km are greater than 10® -.cm. The particle- in-a- 
box formulation is not applicable once the effective size of 

the particle is greater than its distance of random motion. 

Why are random motions damped so effectively in this solid 
rock case? We might expect just the opposite; . 
that the higher coefficients of restitution would give less 
damping. Indeed the damping is very slow: More than 10^ yr 

elapse before relative velocities decrease enough for 
accretion to begin. But precisely because of the high 
coefficient of restitution, accretion can only begin after 
the approach velocities are considerably smaller than the 
1 km bodies' escape velocity. For this reason, the- 
eccentricities and inclinations of the growing bodies are 
small. 

Recognizing that the particle- in-a-box model breaks down 
in Experiment 5 after bodies greater than v40 km have been 
created, we can make reasonable estimates of their subsequent 
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evolution based on our experience with collis'ional .modelling. 
As these '^40 km objects are formed they find themselves 
effectively isolated in nearly circular orbits. They sweep 
up all material which passes within their capture cross- 
section (10® cm) due to Keplerian -differential motion-.. 
Continued growth of smaller bodies left in neighboring zones 
continues to produce more of these 40 km objects. Eventually, 
there are enough of these larger bodies that they begin to 
gravitationally stir one another into more irregular orbits 
(higher e's and' i'.s). Henceforth, the particle- in-a-box 
assumption becomes applicable again. -So long as the random 
velocities do not get too much greater than the larger bodies ' 
escape velocities, accretional growth will then proceed. 
Experiment 6 (Figs. 13, 14 and 15); Weakly Bonded Regolith 

An opposite extreme of material properties was 
introduced by considering parameters appropriate to bodies 
consisting of loosely bonded regolith (Table 1) . Again, the 
population was considered to consist initially of km-sized 
objects. Because of the low value of p compared to previous 
experiments, the initial number of bodies v?as augmented by a 
factor of 4 to keep the surface mass density of the disk 
vlO gm/cm^. The escape velocity of the particles is reduced 
with the density, so we selected lower initial values for e 

and i of 7 X 10“ which is about 3 v . 

. 0 

The subsequent evolution is portrayed in Figs. 13, 14 
and 15. It is quite similar to the general properties of 
previous experiments. The growth of bodies hundreds of km 
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in diameter occurs in about 1000 yr. " This’ rate 'is much faster 
than for rock (Experiment 5) or for our intermediate strength 
material (Experiment 1 ) , presumably because of the increased 
number of initial bodies in this case. It is slov/er than the 
case with initially larger numbers of particles (Fig. 9) . The 
evolution begins with the 1 km bodies colliding, cratering , 
and accreting onto one another. The cratering generates small 
bodies and accretion creates many 4 km bodies within 30 yr. 

Any body smaller than 4 km continues to crater and accrete 
any other body of its own size that it hits. If it hits 
a smaller body, it is cratered and loses some mass and the 
smaller body shatters and escapes. The dominant -process 
is accretion. When two bodies of equal size greater than 4 km 
collide, they shatter but remain gravitationally bound as a 
single object. If a body larger than 4 km hits a smaller 
body, the smaller body is shattered and accreted while the 
larger one is cratered with some ejecta escaping. _ The total 
mass of particles smaller than 30 m, which our program 
neglects, is less than 0.1% of the total. The pattern of 
collision outcomes (Fig. 15) is quite different than for 
previous cases, but the size distribution evolution (Fig. 13) 
follows a similar pattern. Bodies of hundreds of kilometer 
diameter are formed while most of the mass of the system 
remains in 1 km objects. The random velocities for the largest 
bodies are less than for the smaller ones , but all are of the 
same order as the escape velocity corresponding to diameter 1 km. 
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IV. DISCUSSION 


The results just described are for only a few cases from 
the range of possible initial conditions and model parameters 
that one might wish to study. Moreover we have not yet incor- 
porated some physical processes that we expect will be important, 
at least for some cases . Yet the results demonstrate dramatically 
the efficacy of planetesimal accretion. In this section, we discuss 
the significant consequences of our results and some areas in 
which we are continuing the work. 

While there have been uncertainties about many stages of 
solar system origin and planetary accretion, one of the most 
serious questions has concerned the intermediate phase of accretion, 
•i.e.,, growth after the hypothetical formation of . planetesimals from 
gravitational instability (Safronov 1972; Goldreich and Ward 
1973) but before the late stages of accretion when the largest 
bodies have substantial gravitational cross-sections (cf. 

Hartmann and Davis 1975) . It has not been clear how planetesimals 
could have efficiently accreted one another. Our modelling, 
based on detailed physical experiments involving low-velocity 
collisions among rocky bodies, demonstrates that accretion 
through this intermediate size-range is efficient and rapid. 

It is a natural result of low-velocity rebound phenomena dis- 
cussed by Hartmann (197 8). 

Some additional physical processes that might be important in 
this- intermediate stage have not yet been incorporated into our 
numerical simulation. For instance, resonant phenomena might 
accelerate or retard growth in certain zones . Another influence 
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of interest is gas-drag, which seems even more relevant given the 
short time scales (10**' years) in which we are getting substantial 
growth with our present algorithm. Possibly the influence of 
gas-drag on accreted bodies would not differ greatly from the 
influence of the swarms of small planetesimals remaining at 
the end of our simulations, but we intend to model gas. explicitly 
in the future. A potentially disrupting influence on accretion 
would be high-velocity bodies, perhaps scattered into the zone 
of interest by an early-formed Jupiter. 

Later Stages of Planet Growth 

It is interesting to speculate on how later stages of 
planetary accretion might proceed, given the size- and velocity- 
distributions at the end of our simulations. Note that despite 
the development of 500 to 1000 km diameter objects, the bulk of 
the mass in the system remains in the 1 to 2 km diameter bins. 
Derivation of a similar result has been attributed to Y. Nakagawa 
by Hayashi ^ (1977) . This state is similar to 

distributions used implicitly by several workers (Hartmann 
and Davis, 1977; Hayashi et ^. , 1977; Weidenschilling, 1974) . 
as starting conditions for modelling the final stages of 
planetary growth ; a few seed planets with most of the 
mass in a cloud of much smaller particles. Alternatively, 
before seeing our results , one might have imagined 
the intermediate stages to have been characterized by 
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such rapid growth of the smaller bodies that the largest 
bodies grew in numbers fast compared to their growth 
in size. If that were the case, the bulk of the mass 
would have resided in the larger bodies and the later 
stages of planet growth v?ould have involved their mutual 
accretion, rather than accretion of planetesimals by 
seed planets .hundreds of km in diameter. (A size distribu- 
tion with most of its mass in the larger bodies is 
observed today in the asteroid belt, but this is probably 
a product of comminution rather than accretional evolution.) 

^hile it is plausible that the first- formed multi- 
hundred km bodies will act as seeds for subsequent growth 
of full-sized planets, our present model cannot follow 
the detailed processes of such continued collisional 
evolution. This is because the random motions of the 
largest formed bodies have become, at this stage, smaller 
in extent than the dimensions of their gravitationally 
enhanced collision cross-sections. Thus our basic 
particle- in-a-box model breaks dov7n because these bodies 
do not sweep through a representative sample of the entire 
population at a speed governed by random velocities. Instead, 
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a large body is encountered by, and accretes, only those 
smaller objects remaining with orbital semi-major axes close 
enough to the large body's that .differential Keplerian 
velocities bring them within the collision cross-section. 

As this zone is cleared by' accretion, it is conceivable. 
that the large body ' s cross-section would grow so as to 
dominate an increasingly large annulus of planetesimals . 

Alternatively, the large body may become effectively 
collisionally isolated from the rest of the system due to its 
nearly circular orbit. Such isolation would be only temporary. 
Diffusion by interactions of small-scale planetesimals from 

> 

adjacent zones might tend to feed material into the large body's 
accretion zone. If that mechanism is slow or ineffective, 
continued collisional evolution among planetesimals in other 
zones would grow other. 500 km scale bodies by the same process 
which led to the isolated first generation of large bodies. 
Eventually, there would be enough of these large bodies that they 
would begin to perturb one another onto more eccentric orbits 
providing access to one another and to any remaining planetesimals. 
Since the relative velocities due to stirring would be of the 
order of the large bodies' escape velocity, collisions thus pro- 
moted would probably result in accretion (Hartmann 1978) . The 
gravitational scattering of planetesimals from the region of 
the first-formed full-size planet would also tend to break any 
isolation of 500 km objects-. 
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Influence of Orbital Resonances 


In addition, orbital resonances with the first-formed planet 
would have acted to break any such isolation by preferentially 
enhancing the orbital eccentricity of the larger bodies at cer- 
tain semi -major axes, as consideration of the nature of these 
resonances will show. An orbital resonance occurs when a 
particle's orbital period is near a small whole-nxnriber com- 
mensurability with the period of the perturbing planet. 

Repetitive mutual configurations induce a forced eccentricity 
in the particle's orbit, the magnitude of which increases 
with decreasing distance from the exact commensur ability 
{cf. Greenberg 1977) . Similar effects ("secular resonances") 
occur near commensurabilities of precession periods. 

The theory of resonances is a well-studied area of 
celestial mechanics so that computation of forced eccentricity 
is a straightforward procedure. However, a significant forced 
eccentricity at a given semi-major axis does not in itself 
imply enhanced relative velocities, because close particles 
undergo coherent perturbations: the apsides corresponding to 

the forced eccentricities are aligned in such a way as to mini- 
mize collisions. Those particles in a very narrow band near the 
exact resonance have large enough forced eccentricity that their 
radial excursion reaches particles whose motion is not coherent 
with their own. The particles in this narrow band transfer 
random motion to other particles in the vicinity through colli- 
sions. Such collisions could rapidly deplete the population 
of resonant particles unless new material is fed into the 
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resonance zone. This material might be either the scattered 
products of the collisions or material which has undergone 
secular variation of semi -major axis by drag or radiation effects. 

The coherence of forced radial oscillation also breaks down 
at the sudden phase transition across the semi-major axis which 
corresponds to an exact resonance (Greenberg 197 8) . But this 
effect, too, involves only those particles extremely close to 
the critical semi-major axis. 

On the other hand, different sized particles in a popu- 
lation do not have coherent resonant oscillation because, 
while the larger bodies ' orbits respond to resonances and 
achieve an appropriate forced eccentricity, smaller bodies’ 
orbits are drastically, discontinuous ly, and frequently modi- 
fied by collisions with and close approaches to bodies of 
comparable or greater size. Hence, the smaller bodies cannot 
respond to long-term resonant perturbations. In this^way 
forced eccentricities introduce a relative velocity between 
particles in different size regimes. 

The distinction between the response of small and large 
bodies to resonant perturbations can be compared to the 
response of a mass hanging from a spring to a small periodic 
force close to its natural frequency. The resonant amplitude 
can be achieved only if the driving force operates for many 
periods and not if the position and velocity of the mass are 
frequently and arbitrarily re- initialized. These cases would 
be analogous to behavior of the larger and smaller bodies , 
respectively. Note that the larger bodies have their radial 
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oscillations damped by drag due to collisions with the small 
ones. The result is a phase shift and amplitude limits just 
as would occur if drag were introduced to the mass-on-a-spring 
analog. These ideas are explored in more detail by Greenberg 
(1978) . . 

Resonances are thus seen to provide an important mechanism for 
breaking the isolation of larger bodies during the accretion pro- 
cess due to- their nearly circular 'orbits . On the other hand, 
the high relative velocities might have led to catastrophic 
fragmentation rather than accretion at these positions. Perhaps 
growth was favored just adjacent to the resonance positions 
where collisions were reasonably frequent but velocities were 
not too high. . ' 

Several properties of the present planetary distribution 
suggest that an accretional model governed by resonances may 
be relevant. The asteroid belt spans orbital radii which 
correspond to the important low-order commensurabilities with 
Jupiter's orbital period; planetesimals in the belt never grew 
to diameters much greater than 1000 km. j]Chapman and Davis (1975) 
argue that, had they ever exceeded 1000 km, they would still 
survive.] The density distribution within the belt appears 
to be governed by resonances, with either gaps or 
concentrations at commensurable distances. In the outer 
solar system there are striking near-commensurabilities' between 
adjacent planets (Wilkins and Sinclair 1974) ; satellite systems 
contain a statistically significant excess of resonances 
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(Goldreich 1965); and the structure of Saturn's rings appears 
governed by resonances with other satellites (Franklin and 
Colombo 1970) . The terrestrial planets do not exhibit such 
striking mutual commensurabilities , but this might be ex- 
plained by the shift in resonance positions which would have 
occurred in the presence of the early inner disk -of material (more 
dense than in the outer solar system) , just as resonances may 
be shifted in Saturn's rings according to the theory of 
Franklin and Colombo (1970) . 

Earliest Growth 

Although we have applied our model to the intermediate 
stage of planet growth, it may also be relevant for earlier 
stages. In one test case, we applied our model to a case of 
mutual interaction in a population initially of all 1 cm 
particles. Bodies approaching 30 meters in diameter formed 
in only a few years. As in most of our other n\imerical ex- 
periments, most of the mass of the system remained in the 
initial-size particles at the time our particle- in-a-box approach 
became invalid, so it remains to be seen whether direct particle ' 
collisional interaction might be competitive in timescale with 
gravitational instability mechanisms for forming km-scale 
planetesimals . 

Astrophysic'al Applications 

Many astrophysicists have supposed that planet formation 
is a common process in the universe given the dusty clouds 
observed around newly formed stars. Yet so long as there 
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have remained obstacles to modelling the accretion of dust 
into full-sized planets, there has remained the possibility 
that the sun's planetary system is the result of unusual 
circumstances and that other planetary systems are rare. Thus our 
success in attaining rapid accretion through the difficult inter- 
mediate size ranges increases somewhat our expectation ' that 
planetary systems formed around some other stars. 

While the numerical results reported here have concentrated 
on plausible early solar system models, we are currently 
broadening the range of input parameters to discover what 
conditions limit planet growth in the general stellar case. 

For example, one set of numerical experiments demonstrates 
conceptually how a sufficiently high velocity regime may completely 
inhibit growth of a planetary system and produce only a swarm 
of asteroid- like bodies- The experiments indicate that rock 
fragmentation will produce debris extending in size down to 
the size of the homogeneous grains in the shattered material. 

Our work suggests that the collisional evolution of planetesimals 
might produce abundant ym-scale particles which could be driven 
into interstellar space by radiation pressure (Soter, Burns, 
and Lamy, 1976 ) . Thus planetesimal systems would be sources 
of observed interstellar grains as earlier suggested by 
Herbig ( 1970 ) and Hartmann ( 1970 ) . Further details of this 
work will be reported in a future publication in preparation. 
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V. COMCLUSIOW 


Our simulation is still undoubtedly a long way from complete 

reproduction of the collisional evolution of the early solar 
system. The list of physical mechanisms not incorporated 

in the model is presumably endless, although we must assume 
that most would have negligible effect on the results . Our 
algorithm may serve as the basis for testing the degree of 
importance of various phenomena. Certainly, the effects 
of gas drag, orbital resonances, and material scattered 
by Jupiter must be considered. As discussed in the intro- 
duction, our coupled treatment of evolution of mass and 
velocity distributions was largely motivated by the need to 
incorporate these phenomena. Our program is thus structured 
to permit such incorporation- The program is also designed 
to permit updating and refinement of the treatment of impact 
phenomena as more theoretical and experimental work is 
done. Inclusion of some other potentially important 
properties of the population may require structural modifi- 
cation to the algorithm. For example, surface regoliths 
and body rotation rates would evolve synergistically with 
size and velocity distribution during the collisional phase 
of planet formation (cf. Hartmann, 1978, regarding the 
relation with regoliths; Harris, 1977, regarding rotation), and 
these processes will be incorporated • in our program in the future. 

Besides providing a basis for future investigation of 
the relative importance of various phenomena, the simulation 


- 46 - 



is already a higher order approximation of collisional 
evolution than any constructed before. The main conclusions con- 
cerning growth of planets are the following: (a) Collisions beginning 

with km sized bodies rapidly produce a substantial number 
of 500 to 1000 km bodies. This result is based on an experi- 
mentally motivated model of impact outcomes . It requires no 
ad hoc sticking mechanism. (b) The bulk of the mass remains 
as km sized bodies even when the larger objects are formed. 

(c) Random velocities are of the order of the escape velocity 
of the original bodies, not of the large bodies. This result 
contrasts strikingly with the often-quoted conclusion of 
Safronov that velocities were on the order of the largest 
bodies’ escape velocities. Safronov’s result depended on the 
assvimption of (i) a ,power-law size distribution with most 
mass in larger bodies and (ii) an equilibrium velocity 
solution. Our model is independent of such assumptions. In 
fact, neither assumption is justified since we show 
that most of the mass remains in small bodies 
and the growth of large ones occurs too quickly for 
equilibrium to be achieved. (dj Random motion is less for 
big bodies than small ones, because the big bodies 
experience drag due to the smaller ones, while the small 
ones are stirred by gravitational and collisional inter- 
actions with one another. (e) The creation of a small 
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number of bodies in excess of 500 km in a swarm still 
dominated in numbers and mass by much smaller objects 
suggests a picture of subsequent evolution in which the large 
bodies form seeds for the final stages of accretional growth. 
P.urther modelling which properly follows the statistical 
behavior of this small number of large bodies in terms of 
accretion and mutual interaction is needed to continue the 
study of collisional evolution through the formation of full- 
sized planets. 


ACKNOWLEDGMENTS 

This work was supported by NASA Grant NSG 7201 and by 
NASA Contract NASW-2909 with the Planetary Science Institute. 
Conversations with Dr. Donald R. Davis helped develop 
several important ideas. Computation was performed with 
the University of Arizona's CDC-6400 computer. This is 
P.S.I. Contribution Nvimber 83. 


- 48 - 



REFERENCES 


Chapman, C.R. and Davis, D.R, (1975) . Asteroid collisional 

evolution; evidence for a much larger earlier population. 
Science 190 , 553-556. 

Franklin, F. and Colombo, G. (1970) . A dynamical model for 
the radial structure of Saturn's rings. ' Icarus 12 , 

338-347. 

Fujiwara, A., Kamimoto, G. and Tsukamoto, A. (1977). Destruction 
of basaltic bodies by high velocity impact. Icarus 31 , 
277-288. 

Gault, D.E., Shoemaker, E.M. and Moore, H.J. (1963). Spray 
ejected from the lunar surface by meteoroid impact. 

NASA Tech. Note D-1767 . 

Goldreich, P. (1965). An explanation of the frequent occurrence 
of commensurable mean motions in the solar system. Mon. 

Not. Roy. Astron. Soc.' 130 , 159-181. 

Goldreich, P. and Ward, W.R. (1973) . The formation of planetes- 
imals. Astrophys. J. 183 , 1051-1061. 

Greenberg, R. (1977). Orbit-orbit resonances in the solar 
system: varieties and similarities. Vistas in 

Astronomy 21 , 209-239. 

Greenberg, R, (1978). Orbital resonance in a dissipative 
medium. Icarus , in press. 


- 49 - 



Greenberg, R. , Davis, D.R., Hartmann, W.K. and Chapman, C.R. 
(1977) . Size distribution of particles in a planetary 
ring. Icarus 30 , 769-779. 

Harris, A.W. (1975) . Collisional breakup in a planetary ring. 
Icarus 190-192. 

Hartmann, W.K. (1969) . Terrestrial, lunar and interplanetary 
rock fragmentation. Icarus 10 , 201-213. 

Hartmann, W.K. (1970) . Growth of planetesimals in nebulae 
surrounding young stars. In Evolution Stellaire avant 
la Sequence Principale 59 , pp. 215-227, Univ. Liege Collog. 

Hartmann, W.K. (1978) . Planet formation: mechanism of 

early growth, Icarus , in press. 

Hartmann, W.K. and Davis, D.R. (1975), Satellite-sized 
planetesimals and lunar origin. Icarus 24 , 505-515. 

Hayashi, C., Nakazawa, K. and Adachi, I. (1977).. Long term 
behaviors of planetesimals and the formation of the 
planets. Kyoto Univ. preprint KUNS 379. 

Herbig, G.H. (1970). Introductory remarks. In Evolution 
Stellaire avant la Sequence Principale 59 , pp. 13-26, 

Univ. Liege Colloq. 

Isaacman, R, and Sagan, C. (1977) . Computer simulations of 
planetary accretion dynamics: sensitivity to initial 

conditions. Icarus 31 , 510-533. 

Kaula, W.M. and Bigeleisen, P.E. (1975). Early scattering 
by Jupiter and its collision effects in the terrestrial 
zone. Icarus 25, 18-33.’ 


- 50 - 



Kuiper , G.P. (1951a). On the origin of the solar system. 

In Astrophysics (J.A. Hynek, ed) , Chapter 8, McGraw-Hill, 

New York. 

Kuiper, G.P. (1951b). On the origin of the solar system. 

Proc. Nat. Acad. Sci. 37 , 1. 

Kuiper, G.P. (1974) . On the origin of the solar system. 

Celestial Mechanics 9_, 321-348. 

Lecar, M. and Franklin, F.A. (1973) . On the original distribu- 
tion of the asteroids. Icarus 20 , 422-436. 

Marcus, A.J. (1969) . Speculations on mass loss by meteoroid 
impact and formation of planets. Icarus 11 , 76-87. 

Moore, H.J. and Gault, D.E. ,(1965) . Fragmentation of spheres 
by projectile impact. Astrogeologic Studies , U.S.G.S. 

Annual Progress Report, 127-150. 

Safronov, V.S. (1972). Evolution of protoplanetary clouds 
and formation of the earth and planets. NASA TT F-677 . 

Soter, S., Burns, J.A. and Lamy, P.L. (1976). Radiation 

pressure and Poynting-Robertson drag for small spherical 
particles. Proc. Internatl. Astr. Union Colloq. 39 . 

Stoffler, D., Gault, D.E., Wedekind, J. and Polkowski, G. 

(1975) . Experimental hypervelocity impact into quartz 
sand: distribution and. shock metamorphism of ejecta. 

4062-4077. 

Ward, W.R. (1976) . The formation of the solar system. In Frontiers 

of Astrophysics , (E.H. Avrett, Ed.), Harvard Univ. Press, pp. 1-39. 

Weidenschilling, S.J. (1974). A model for accretion of the 
terrestrial planets. Icarus 22 , 426-435. 


-51 



Weidenschilling , S.J. (1975). Mass loss from the region of 
Mars and the asteroid belt. Icarus 26 , 361-366. 
Wilkins, G.A. and Sinclair, A.T. (1974). The dynamics of, 
the planets and their satellites. Proc', R. Soc. 

London A. 336, 85-104. 


- 52 - 



TABLE 1; 


PARAMETERS FOR VARIOUS EXPERIMENTS, 
ALL UNITS ARE CGS , 


Experiment 


# 1 , 2 , 3 , and 4 #5 


Comments ; 


Intermediate 

Material 


Solid 

Rock 


Max. Vel. for 100 4000 

Rebound w/o Cratering 


Coeff. of Restitution 0.5 0.86 

r, c . 

1 


Modified Coeff. of Rest. 0.1 

c . 

ii 

Mass Excavated Coeff.' 10~® 

K 


Ejecta Vel. Coeff. 3 x 10^* 

C . {cf. Fig. 1) 

S3 

Density 3 

P 


Impact Strength 
S 


3 X 10 ^ 


0.86 


10“9 


2 X 10 ^ 


4 

3 X 10 ^ 


#6 

Loosely Bonded 
Regolith 

. 1 

0.01 
0.001 
3 X 10 ~^ 

10 ^^ 

0.9 

lO'* 


*NOTE; In #2 and #3 all crater ejecta v?as given a velocity of 
0.005 times impact velocity. 



FIGURE CAPTIONS 


Figure 1: The ejecta velocity distribution of Gault ^ al . 

(1963) for basalt (heavy curve) is approximated by a straight 
(dashed) line. Results of Stoffler et (1971) suggest 

velocities are a factor of 100 smaller for sand. The dotted 
segment is our extrapolation of the sand line beyond the range 
measured by Stoffler et a]^. We consider an intermediate distribu- 
tion, with coefficient c . = 3 x 10® cgs, as well. 

Figure 2: Mapping of outcomes as function of impact velocity 

and mass ratio. Velocity is upper limit for rebound. 
Catastrophic fragmentation (shattering) occurs if impact 
strength (critical energy/unit volume) is exceeded. 

Figure 3: Schematic representation of our modeling of 

change in mass of an impacted body a's a function of impact 
velocity and energy- delivered to the body . Compare with 
Hartmann’s experimental results (Hartmann, 1977, Figure 1). 

Figure 4: Particle size distribution as a function of time 

for a material intermediate between loosely bonded regolith 
and solid rock. Experiment #1. 

Figure 5; Matrix indicating results of collisions between 
pairs of bodies of various sizes in Experiment #1. 

Code: 1 = Escape after outcome (i) ; 

2 = Both bodies cratered (ii) and escaped from one another; 

3 = One body shattered, its debris escapes other body; 

4 = Both bodies shattered, debris escapes; 

5 = Accretion after outcome (.i) or (iv) ; 



6 = Accretion after outcome (ii) ; 

7 = Accretion after outcome (iii) . 

Figure 6: Eccentricity distribution as a .function of time and 

particle size for Experiment #1. Inclinations are similar. 

Initial value, e = 5 x 10“**, is shown by +. Corresponding 
random velocities are shown on right hand scales in terms 
of 1 km escape velocity, 20 cm/sec. 

Figure 7: Size distribution evolution for Experiment #2 with 

crater ejecta escape effectively suppressed but otherwise 
parameters and initial conditions are the same as for Experiment 
#1. Note similarity of growth. 

Figure 8; Size distribution for Experiment #3 which was identical 
to Experiment #2, but with finer size resolution. Results are 
similar indicating that they are not limited by our coarse 
size bins . 

Figure 9: Size distribution evolution for Experiment #4 which 

was identical to Experiment #1, except with 100 times as many 
initial bodies. Principal difference is contraction of the 
time scale. 

Figure 10: Size distribution for a case (Experiment #5) using 

impact parameters appropriate for solid rock. 

Figure 11: Collision outcome matrix for Experiment #5. See 

caption of Fig. 5 for code. 



Figure 12: Velocity distribution for Experiment #5. Here, 

V 25 cm/sec ; 
e 

Figure 13; Size distribution for a case (Experiment #6) , 
using impact parameters appropriate tO’ weakly, bonded regolith. 

Figure 14: Collision outcome matrix for Experiment #6.. See 

caption of Fig. 5 for. code. 

Figure 15: Velocity distribution for Experiment #6. Here 

V 'v 10 cm/sec. 
e 
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